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We present results on a stochastic forest fire model, where the influence of the neighbour 
trees is treated in a more realistic way than usual and the definition of neighbourhood 
can be tuned by an additional parameter. 

This model exhibits a surprisingly sharp phase transition which can be shifted by 
redefinition of neighbourhood. The results can also be interpreted in terms of disease- 
spreading and are quite unsettling from the epidemologist's point of view, since variation 
of one crucial parameter only by a few percent can result in the change from endemic to 
epidemic behaviour. 
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1. Introduction 

Cellular Automata (CA) are recognized as an important modelling tool in many 
fields of science, including also biology and especially ecology. The main character- 
istics of the CA approach are summarized in three rules (see [T]): 

(Rl) Space and time are discrete: There is a (often one- or two-dimensional) grid of 
cells that is only viewed at distinct (equidistant) timesteps. 

(R2) At each timestep, each cell has one (and only one) state taken from a limited 
number of possible states. 

(R3) There are simple but universal deterministic update rules: The state of a cell 
at a given timestep depends only on its own state and the cell states in its 
neighbourhood, all taken at the previous step. These rules are the same for all 
cells and timesteps. 
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Fig. 1. Different neighbourhoods for Cellular Automata: (a) von Neumann neighbourhood, (b) 
Moore neighbourhood, (c) extended Moore neighbourhood and (d) neighbourhood mediated by 
the Moore parameter ttm 



A more detailed discussion of the CA approach and its history can be found for 
example in [2J. CA models which strictly respect these three rules (and, in many 
cases, make use of only two cell states) have been intensively studied and sometimes 
can give considerable insight in the foundations of complex behaviour. 

Each rule can be modified in order to extend the range of systems which can 
be simulated. We soften rule (i?3) by making the update rules still universal, but 
stochastic, employing (pseudo) random numbers to determine cell evolution. 

A topic that cellular models often excel in is ecological modelling, especially the 
study of disease spreading. In fact, there is a wide variety of systems that can be 
modeled with CA, and disease spreading models of this type have been successfully 
applied to plants (e.g. [3]), animals (e.g. [4]) and, in part, human beings (e.g. [5]). 

Often CA models are not applied to special (realistic) cases, but general prop- 
erties are studied, as for example in [3] or [TJ. One of the most interesting (and 
most dramatic) features of CA models is the possibility of "phase transition" (see 
for example [5J[H])- 

In the following we will study such transitions as well - where the term phase 
transition is used in a more general meaning than the one common in physics since 
there does not necessarily exist any quantity analogous to free energy. Still, some 
order parameter can be found which allows to clearly distinguish distinct phases of 
the system and on the border there are indications for critical behaviour. 



2. The Definition of Neighbourhood 

On a 2D lattice, the neighbourhood of a cell may be defined in many ways; the most 
popular ones are, however: 

• von Neumann neighbourhood: The neighbourhood consists only of the four ad- 
jacent cells. 

• Moore neighbourhood: In addition to the von Neumann neighbourhood, also the 
four next-nearest-neighbours are taken into account. 

• extended Moore neighbourhood: An even larger number of cells is considered as 
neigbourhood. 
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These types of neighbourhood are illustrated in figure [1] The neighbourhood type 
and the number of cell states determine the number of rules one has to specify (in 
principle). For a neighbourhood of N cells and s possible cell states one needs s N+1 
rules. (This means, for example, 2 9 = 512 for a two-state model on a 2D grid with 
Moore neighbourhood as for example Conway's Game of Life.) 

Therefore in many models, the updates rules depend only on the number of 
neigbours in a certain state, but not on their actual position. This of course greatly 
simplifies the speficication of update rules. 

When doing this, an interesting type of rules are those that use an extended 
neighbourhood (at least Moore-type), but reduce the importance of non-adjacent 
cells for the update. For update rules that only depend on the number of neigbours 
in a certain state, their number N can be redefined as 

N := iVadj + 7TM • iVnoa-adj (1) 

where iV a dj denotes the number of adjacent and iV non _ a dj the number of non- 
adjacent cells that are nevertheless regarded as neighbours. In case of the next- 
nearest- neighbours, the Moore parameter ttm mediates the transition from a pure 
von Neumann neighbourhood (ttm = 0) to a full Moore neighbourhood (ttm = 1). 
The accordingly modified neighbourhood is also sketched in figure [TJ Models of this 
kind have been studied in this paper, and the influence of ttm on certain results has 
been checked. 

3. The Forest Fire Model 

Our model is a variation of wide-spread forest fire models (see p~0 ^ ITT1 . fl2| . fT3| . IT4 l fl5 | 
116] , which can also be re-interpreted as a disease-spreading model. In this article, 
however, the vocabulary will be mostly the one of living, burning and burnt-out 
trees. 

3.1. Definition of the Model 

The forest fire simulation takes place on a quadratic two-dimensional grid (in all 
the following simulations 128 x 128) with periodic boundary conditions and a medi- 
ated Moore neighbourhood as described in section [21 The states were encoded the 
following way: 



X 


colour 


state of tree 


-1 


black 


dead (burnt-out) tree 





white 


empty cell (or fire-resistant tree) 


1 


green 


living tree 


2 


red 


burning tree 



The only possible transitions are those from living to burning state, [1 — > 2], and 
those from the burning state to the burnt-out one, [2 — > (— 1)]. The probability that 
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Fig. 2. (a) Sketch of the basic forest fire model (b) Expamples for weak, drastic and critical 
decreases of the number of burning trees in a smoothed time series 



a living tree at (i, j) catches fire is given by: 

P[i^2] (i, j, 1 1 M-i) = i ■ {1 - (1 -Pbf^- 1 } 



(2) 



where pb denotes a fixed burning probability. The state of the cell's neighbourhood 
at the previous timestep is symbolized by the condition complex A/t-i, it explicitly 
enters the formula as Nb,t-i, the number of burning neighbours at (t — 1), 



N, 



b,t-l 



(3) 



A^e-f-iAi} Ai,je{-i,o,i} 

|Ai|+|A 3 -|=l |A i |+|A ;f |=2 

A burning tree is converted to a burnt-out one after one timestep, formally 

P[2->-l] («, j, *|M-i) = ^x i>3 ,2 

A symbolic sketch of the model is given in figure [U a. Equation @ together with 
the definition of Nb t t-i expresses the increased probability of a tree catching fire 
if it has more than one burning neighbour. In the case N^t-i = 1 the probability 
reduces to 



(4) 



P[i~>2] {xi,j,t\J\f t -i) = S Xitjt i •{!-(! ~Pb)} = - Pb 



(5) 



3.2. A Typical Simulation Run 

At the beginning of the simulation, living trees are placed with a density d s tart G 
[0, 1] while the remaining cells stay empty. Than a few number of burning trees 
are inserted - starting with only one burning tree would dramatically increase the 
probability of a " false start" . 
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Fig. 3. Snapshots of a forest fire simulation and corresponding time series of burning trees (critical 
case, p b = dstart = 0.63) 

An example of a typical simulation run is shown in figure [3] for the choice pt = 
^start = 0.63. Even some runs with varied parameters already indicate that there 
is a relatively sharp "phase transition" between those configurations where the fire 
dies out within a few steps and those where almost the whole forest is affected. 

This is shown in figure 0] where the end configurations of different characteris- 
tic simulations runs are given. There is indeed an analogy to phase transitions in 
physical systems - for example the Ising (see figure [SJ taken from p~7] ) . 

In the Ising model, for temperatures below T c w 2.268254878 the system is quite 
homogeneous (ordered with only small fluctuations). It is quite homogeneous as 
well for temperatures significantly above T c (disordered with only small correlated 
regions). At T = T c , however, the picture is drastically different: There are clusters 
on all length scales and the correlation length diverges. (For all simulations on a 
finite grid, the correlation length cannot diverge; thus the picture of clusters at all 
length scales remains also true for T w T c .) 
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(a) (b) (c) 




Fig. 4. Typical end configuration and time series of burning trees in forest fire simulations: (a) 
subcritical case (p;, = d s tart = 0.60), (b) critical case (pj, = dstart = 0.63) and (c) supercritical 
case (p;, = <4tart = 0.66). Note the different timescales, in particular the fact that in case (b) the 
fire lasts about twice as long as in case (c). 




Analogously we denote cases where the fire dies out within a few timesteps and 
leaves behind a homogenous, almost unaffected system as subcritical, as illustrated 
in [Ha. A situation as depicted in figure [Hb where one finds burnt and unaffected 
clusters at all length scales is classified as critical. If, as depicted in figureHJc, the fire 
affects the whole system and leaves behind again a homogeneous system of burnt 
tree we denote this as supercritical. 
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3.3. Relation to Other Models 

In contrast to the majority of the models mentioned in the introduction to sec- 
tion [3J our model does not exhibit self-organized criticality. For example in the 
Clar-Drossel-Schwabl-model [12] the combination of a growth process and repeated 
ignition of forest fires is believed to push the system into a critical state. From 
the exposition given in subsection 13.21 this is easy to understand on an intuitive 
level. As long as the system is in a subcritical state, fires will die out after a few 
timesteps, resulting in net growth. In a supercritical state, large clusters will burn, 
which yields net loss of trees. The combination of both mechanisms will lead to a 
critical quasi-equilibrium. 

In our model no growth mechanism is present, so - as it is the case in typ- 
ical statitical mechanics systems - parameters have to be fine-tuned to access a 
critical state. Accordingly this model provide^ a very close look at the bruning 
process of a cluster, a "magnifying glass" which allows to study certain properties 
of self-organized criticality whithout the need to be very close to the critical point. 
Otherwise the critical state is hard to access numerically, since in order to deal 
with diverging correlation lengths, huge grids are required to have finite-size effects 
under control. 

As outlined in [5] the main reason why this model has been studied is the 
fact that it can easily generalized to a semi-realistic disease-spreading model. A 
publication on that model is in preparation (see also comments in section [6]). 

4. Analysis Tools 

After having performed the simulation (which is finished when there are no burning 
trees left), it is the main task to extract the essential information and to condense 
it to a small set of significant numbers that can be easily compared for various runs 
and different parameter configurations. 

There are two sources of information that are easily available: the final tree 
configuration and the time series of burning trees (both depicted in figure [4|) . The 
main quantities which can be extracted from them are summarized in table [9] 

4.1. Cluster Analysis 

The methods covered in this section are in a way related both to determining 
the fractal dimension of geometric objects by box-counting and the ideas of the 
renormalization group [18) , in particular when formulated in the language of block- 
spin transformations a la Kadanoff. 

Starting point is a the end configuration after a simulation run, i.e. a lattice 
of cell states := {sij} (as they were defined in equation ([3]): s = —1 for a 
burnt-out, s = for an empty and s = 1 for a living cell). Now multiple cells are 

a for burning probability pi, = 1, if follows the path of the most popular models 
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Fig. 6. Combining cells as described in equation J7J 



combined to new one, 

s i, j = f { s iij,l, jij.l i s iij,2,jij,2 > S Uj,3,3ij,3 J ■ ■ • ' S iij,k,jij,k ) ' (6) 

a procedure that yields a smaller lattice of new cells with states = {s^ j}. This 
process can of course be repeated, resulting in a sequence of lattice configurations: 



where is a lattice consistent of too few cells to repeat the procedure defined in 
equation ([6]). Here, four cells were combined to a new one via 

s 'i, j = S ig n ( s 2i-l, 2j-l + s 2i-l,2j + s 2i,2j-l + s 2i,2j) (J) 

a process which is illustrated in figure [6] Of course, other procedures are possible 
as well, and others are necessary to treat models with other types of states. 

A suitably chosen characteristic observable 0(n) can be evaluated for each of these 
configurations, 

0(n) = O(S^), (8) 

and its values may give insight to the typical length scales present in the system. 
On this behalf, one can try to fit an analytic function O(x) to them and read off a 
the size distribution of clusters. The most obvious choice for O is 

0W = j (9) 

(with N n denoting the total number of cells in a configuration S^) which was 
employed here, but alternatives for O are possible as well. 
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Of course, after a sufficient number of recombination steps, in most cases values 
of +1 or —1 will remain. The first case corresponds to the scenario where the fire has 
stopped after a comparatively small number of timesteps (subcritical), the second 
is valid for mostly burnt-out trees (supercritical). In both cases, also in view of well- 
known results on fractal respectively box-counting dimensions, O(n) from equation 
(|9|) can reasonably be approximated by an exponential function: 

f(n) = 1 - e- a ' n (10) 

for dominance of +1, respectively 

fin) = -1 + e~ a ' n (11) 

in the case of —1 dominating. The constant a e here is a measure of cluster size: 
When a e is small, this means that there are also larger clusters of opposite cell 
state; a large value of a e indicates rapid "saturation" and therefore the existence 
only of small clusters. 

In the critical case, however, living and burnt-out trees are equally present in 
the whole simulation region, so this results in values of zero repeatedly emerging on 
all length scales during the combination process, and so the exponential fit often 
fails 

So in addition, also a linear (least-square) fit has been performed. Of course, in 
most cases, a linear function is only a rather poor description - but it is reliable in 
the sense that there is a straightforward way to find the fit and one always has a 
result (more or less meaningful). So in addition to a e from equation 1 1 01 respectively 
111! also a parameter a\ for a fit 

f(n) = ai n + C (12) 

has been retrieved. Both parameters usually show the same qualitative behaviour. 
(An important exception for pi, — > 1 is discussed in section l5Tl ) 

For the critical case the cluster size distribution is expected to obeye a power- 
law; while this issue has been left to further investigation (which has to be performed 
on much larger grids and to be compared with results for SOC models mentioned 
in section [3]), a sign for an arising power-law ist the fact that in the critical case an 
exponential fit (flOl [TT|) becomes unstable and often breaks down completely. 

4.2. Time Series Analysis 

Information can also be retrieved from the time series of burning trees. Examples 
for such a time series are given in figures [31 HI El an d EJ There is number of mean- 
ingful and interesting quantities that can be extracted from such a time series by 
rather simple means. Some other parts of the analysis done, however, rely on more 
sophisticated techniques, namely discrete Fourier transform (see section T4. 31) . 

b In the following calculations, the fit has been done with the MATLAB nlinfit routine with stan- 
dard parameters (version 7.1.0.183), which reports failure if the required precision cannot be 
achieved. 
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Fig. 7. Time series and cluster analysis for typical runs of forest fire simulations (N = 128): (a) 
subcritical case (pi, = dstart = 0.61), (b) supercritical case (pi, = ritart = 0.66). In each case, the 
results of exponential and linear fit to cluster compression (equations HOI and lll I . the time series of 
burning trees, a corresponding histogram, the Fourier transform, a smoothed version of the time 
series and at last a very rough histogram are shown. 
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(a) 
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Fig. 8. Time series and cluster analysis for two critical forest fire simulation runs (TV = 128, 
Pb = 4tart = 0.63). As in figure [T] exponential and linear cluster fit, the time series of burning 
trees, the Fourier transform, a smoothed time series and two histograms are given. In (b), however, 
the exponential fit (as a tool for cluster analysis) fails. 
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Some characteristic numbers can be read off from the time series without prob- 
lems: the total burning time tburo (timesteps from the first tree catching fire to the 
last one becoming burnt-out) and the maximum number of burning trees iVb U rn,max- 
The three main cases (subcritical, critical and supercritical) are shown in figure [4] 
(where one should mind the different scales, both for time and number of burning 
trees). The corresponding time series is in all cases quite characteristic: 

• In the subcritical case, the fire only affects a small fraction of the system. In 
the strongly subcritical region, it is extinguished after a few timesteps, but for 
parameters approaching the critical region, it may last much longer. However, 
at any timestep the number of burning trees is comparatively small. There is 
of course at least one maximum, but there may be several of them, in any way 
they tend to be not significant. 

• In the critical case, the progression of the fire is far more complicated than in 
the other cases. Large unaffected clusters may form (and sometimes be partially 
consumed later by the fire). The maximum number of burning trees is smaller 
than in the supercritical case (but significantly higher than the subcritical val- 
ues). Due to this and to the emerging cluster structure, the fire tends to last 
longer (up to approximately the double supercritical burning time). 

• In the supercritical case, all but a few trees are destroyed by the fire. Since the 
fire propagates nearly with maximum speed, it dies off soon after the whole 
system has been affected. There may be some small clusters, separated from 
the bulk, that allow the fire to survive some steps longer, but nevertheless there 
is usually a significant maximum and in virtually any case a rapid decrease (as, 
for example, in figure[7|). In the supercritical limit (pb — > 1, d s tart - ► 1), the time 
series approaches a sawtooth with a clear maximu and a sudden decrease to 
zero one or two steps later. 

Also the fraction of burnt trees /burn can easily be calculated. In this special case, 
it is simply the sum over the whole time series, divided by the initial number of 
trees. Again, in the subcritical case, /bum will be small, in the supercritical it will 
approach one, and for the critical case there will be a wide range of possible values 
- even for the same set of parameters. 

Possibly tbum is the best value to indicate the phase transition from sub- to 
supercritical configuration. When starting subcritical and then slowly increasing pb 
and/or d s t a rt, one will notice a dramatic increase in tbum when reaching the critical 
region and a significant decrease when leaving it towards supercriticality. (This is 
discussed in more detail in section f5. II ) 

Also iVburn,max and /burn are possible indicators for the transition, but since they 
usually monotonically increase on the path from sub- to supercriticality, it may be 
harder to identify the phase transition when watching these quantities. 

c On an TV X iV-grid, a simple calculation shows that there is a maximum of iVbum max = 4(A r — 2) 
at t = ~ 2 for an even N respectively iVhum,nuoE = 4(7V — 1) at t = ~ for an odd grid size N. 
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Another sort of information can be extracted from histograms. Figures [7] and [8] 
also include two histograms each - one with 20 bins and one with only three bins. So 
the second one contains three characteristic numbers ho, hi and hi which give the 
fraction of timesteps when there were comparatively few, medium or many burning 
trees H 

A quantity that may become quite interesting as well, is the number of dramatic 
decreases of the number of burning trees - but due to the rather noisy nature of the 
time series this information is relatively hard to extract from the original series. So 
it is postponed to the next section which deals with Fourier analysis. 

4.3. Fourier Analysis 

From the time series one can deduce the Fourier spectrum as well. This has been 
done with the MATLAB (R) FFT routine (using Fast Fourier Transform). The 
Fourier transform is also displayed in figures [7] and In all cases, low frequencies 
dominate the spectrum (which is already clear from the shape of the time series). 

In addition there tend to be several other peaks that represent higher frequencies 
modulating the signal; for very high frequencies, these peaks drown in noise. So on 
the one hand, the positions and intensities of such peaks (or at least the most 
important one) may characterize some aspects of the time series. 

On the other hand, by cutting off the noisy high-frequency part of the spectrum 
and transforming back, one has a smoothed version of the time series (as also 
shown in figures [7] and [5]) . From this modified time series, for example the number 
of decreases can be read off. 

These decreases have been separated in weak (less than 25%), drastic (25% to 
50%) and critical (more than 50%) ones, and the counts of the latter two have been 
regarded as significant for the fire spreading. Of course there will be at least one 
dramatic decrease (when the fires finally dies off), which is trivial and is therefore 
not counted here. Especially in the case of critical fire propagation, there tend to 
be several more such decreases. (See for example figure [2]b.) 

The appearence of drastic and critical decreases seems to characterize critical 
behaviour quite well, since it is more or less absent in sub- or supercritical sim- 
ulation data. This can also be understood from the characteristics of critical fire 
propagation: On the edge between affecting most cells and dying off soon, there are 
repeatedly situations when only a few burning trees are left and those when a whole 
cluster catches fire. When smoothing out the time series, setting the cutoff point 
usually requires some fine-tuning (and it of course may affect all results). 

For the Fourier analysis presented in section [5] the method described so far has 
been slightly modified: In order to increase the comparability of results, short time 
series have been expanded to a length of 2 8 = 256 steps (by adding zeros after the 
simulation data). So one can establish a more unified scale for the Fourier transform. 

d Of course one has the relation ho + hi +fe = 1; thus only two of these numbers are independent. 
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Q 


meaning 


sec. 


/burn 


fraction of burnt trees (burning rate) 


IP1 


^burn 


total burning time 


EJ 


^burn,max 


maximum number of burning trees 


SJ 


a e 


exponential lit parameter from cluster analysis 




Oil 


linear fit parameter from cluster analysis 


EH 


h 


fraction of steps with "few" 
(e [0, jVbur ^- max ]) burning trees 


EM 


hi 


fraction of steps with "several" 

^ iVburn.max ^ 2iV bu , n , max ^ ]-, urnm g treeS 


EJ 


h 2 


fraction of steps with " many" 






(€ ( bur 3 "'"" x . iV burn>max ]) burning trees 


SJ 


VD 


frequency of dominant Fourier peak 




Id 


intensity of peak at vjj 


1431 


^dras 


number of drastic decreases (25% to 50%) 


1431 


^crit 


number of critical decreases (> 50%) 


1431 



Fig. 9. Summary of the quantities Q eventually used as indicators of the fire propagation (and 
especially the phase transition), "sec." denotes the section where the appropriate quantity is ex- 
plained in detail 

From there, all but the lowest 32 frequencies have been cut from the spectrum (that 
contained initially \ -2 s = 128 frequencies). While this choice is of course arbitrary, 
it yields trustworthy results that mostly coincide with those put forth by the best 
instrument for such types of analysis - human eye and brain. 



5. Results 

Here some fundamental results are presented. In contrast to the few (but neverthe- 
less characteristic) examples presented so far, now the complete parameter space of 
the model has been scanned systematically. 

Since there are only two parameters that strongly influence system behaviour 
(namely pb and d s tart) such a scan is possible with moderate computational effort 
and the results can be displayed in a 2D-plot. 

The grid size has little influence on any result (as long as the lattice is sufficiently 
large, so that the initially burning trees incur no significant bias and subcritical 
propagation affects only a small fraction of the system). So we used a 128 x 128 cell 
grid for all simulations. 

The Moore parameter ttm has greater influence on the result, in section 15.21 
however, we will demonstrate that while there are quantitative changes, the quali- 
tative behaviour is not modified. So all other simulations have been performed with 
km = 0.5. 
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The parameters pb and c? s t a rt have been varied from zero to one in steps of 
Apt — Ad s tart = 0.01, for each configuration there have been 40 simulation runs 
where the quantities given in section [4] have been extracted^] 

Both mean values (indicated by a bar, as for example in tburn) and standard 
deviations (here indicated by a capital delta, A) of these quantities are, for all 
configurations, given in figure [T0l 

5.1. The Phase Transition 

As it has already been mentioned, the model exhibits a (surprisingly sharp) "phase" 
transition between two radically different kinds a behaviour, referred to as subcritical 
and supercritical, as described in section 14.21 On most plots in figure 1101 a. the 
transition line is clearly to be seen. 

This line is approximately a hyperbola with pb ■ d s tart ~ C — const where we have 
the C w 0.4 for the constant^! While the hyperbola is at least no bad approximation 
for the transition line, the true form can not be deduced from the given data - it 
may even have a fractal substructure. Its width, however, does not exceed a few 
percent of change in both parameters (and probably it is still widened due to finite 
size effects and the limited amount of data given.) 

The burning rate /bum increases from below 10% to above 90% on this line. 
On the same line in configuration space, the total burning time iburo reaches its 
maximum value, decreasing both towards the sub- and the supercritical region. The 
burning maximum iVbum.max begins to increase from significantly below 100 to more 
than 200, but this transition is not that sharp and there is a further increase up to 
about iVburn,max = 500 for p b — > 1 and d star t — > 1- 

For the exponential decay parameter a e there seems to be no clear transition 
at all, this impression, however, is created only by the plot scales. Since there are 
both regions with > 25 (for d s tart <C 1 and (1 — pb) <C 1) and < —25 (for 
(1 — dstart) -C 1 and (1 — Pb) <C 1), the line of w that indicates the transition 
is hardly visible. 

Since the range of values is more limited for ai, the phase transition is easier 
to recognize for this fit parameter. The transition again has al ~ 0, but now, such 
values are also reached for pb — v 1, in region where a^T takes its extreme values. 
This, however, is clear from the nature of the different fits. In these regions, the 
exponential decay is very fast, so already from the second step on, the data points 
lie approximately on a straight horizontal line, and therefore we have ai « 0. 

The histogram data (ho to /12) is more difficult to interpret. Also in these plots, 
the transition line is visible, but it is only one of several regions where interesting 
changes occur. Also for frequency and intensity of the Fourier spectrum, the inter- 

e The exponential fit parameter ct e has not always been available, as it has already been discussed 
in subsection |4?T] 

'This value of course depends on Wm> as it can be seen in section [572] In the present case, we have 
C ~ 1 — pc, where pc ~ 0.5928 denotes the percolation limit in two dimensions, see |19| . 
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Fig. 10. Significant system characteristics for pt = : .01 : 1, dstart = : .01 : 1 and 40 iterations: 
(a) mean values, (b) corresponding deviations. The quantities displayed are from left to right (i) 
in the first row burning rate /burn; total burning time tburm maximum number of burning trees 
^bum.max an< l exponential decay constant a e , (ii) in the second row the linear decay constant o; 
and the histogram values h\, h,2, h-j, (iii) in the third row the frequence v^> of the dominant peak, 
its intensity Id, the number n(j ras of drastic and n cr it of critical decreases. 
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pretation is not that easy. For the product pb ■ ci s tart being small (i.e. <C 1), one has 
a low maximum intensity at a low frequency. This can be accounted to a very short 
time series, yielding only a few frequencies. (The following zeros affect the scaling, 
but not the frequencies present in the signal.) 

With pb and/or d s t a rt increasing, first frequency and intensity of the peak tend 
to increase, indicating a signal fluctuating an a small scale. Near the critical line, 
however, there is a significant decrease in mean peak frequency, indicating modula- 
tion of the signal on much larger scale. This is also consistent with the increase of 
the number of drastic or critical decreases - one often has several times when the 
fire has nearly died off, but also periods of chain-reaction-like propagation. 

These results (especially when interpreting the model as one of disease propaga- 
tion) are in a way unsettling: Even an small change of either population density or 
susceptibility may introduce the transition from an isolated (endemic) disease to 
an epidemic or even pandemic one. Especially a decrease in resistance can easily 
be induced by natural catastrophes or starvation - the black death in Europe from 
1347 to 1353 or the Spanish influenza after world war I are likely to be interpreted 
this way. 

On the other hand, they also indicate that even small changes (in population 
density as well as in susceptibility) can reduce an epidemic diseases to a mere 
nuisance. However, one should be cautious when transferring the results of such 
models to systems which they are not intended to be described (as it is human 
society). 

5.2. Influence of the Moore parameter itm 

As defined in[2l the Moore parameter ttm is a way to formally mediate the transition 
from a von Neumann to a Moore neighbourhood. In the previous simulations, ttm — 
0.5 has been employed. Now we study the effects of changing ttm- 

Therefore simulations as described in lS.ll have been performed for ttm — (pure 
von Neumann neighbourhood) and ttm — 1 (full Moore neighbourhood). The results 
are given in figure fTTl 

As it is also the case in figure [TQl there is a clear phase transition (indicated 
especially by /burn, itmrn and a;), but the transition line is shifted - towards a larger 
product pt ■ dstart for ttm — 0, towards a smaller pb ■ d s tart for ttm = 1- (This result is 
fully consistent with what could be expected, since a higher ttm results in a higher 
fire propagation probability.) 

Since the behaviour is qualitatively the same for all choices of ttm, it seems justi- 
fied to perform this type of simulation with just one value for the Moore parameter 
(where we have used ttm = 0.5). 
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Fig. 11. Results of changing the Moore parameter: (a) ttm = 0, (b) ttm = 1| (Pb = : -01 : 1, 
rfstart = : .01 : 1, 10 iterations each). The quantities displayed are the same as in figure [TOl except 
that only mean values, not deviations are shown. 
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5.3. Individual Susceptibilities 

While the model presented here is of course limited in its range of applications, there 
is nevertheless a simple modification that can be implemented with virtually no 
additional effort - the substitution of a global burning probability by an individually 
varying one. 

So we define a matrix P with elements pij that describe the individual burning 
probabilities of the cells (i, j) and replace equation ([2|) by 

P[i^2] (i, j,t\K-i) = S Si , jA ■ {1 - (l-ftj^'*- 1 } (13) 

The distribution of the elements pij can in principle be chosen arbitrarily, but it 
seems reasonable to pick them from a Gaussian distribution centered at p b : 

p id =med{0, G(p b ; a b ), 1} = mcd {0, p b + a b G(0; 1), 1} (14) 

with mcd denoting the median value and G(/i; a) a random variable from a Gaussian 
distribution with mean fi and variance a. For a b — this is exactly the model 
already studied, but for a b ^ (where we use without loss of generality a b > 0) 
new effects may arise. So we studied the various cases; the results for a distribution 
(u b = 0.1) and a broad distribution (a b — 0.5) are displayed in figure \Y2\ and figure 
[T3l (Compared to figures fTUl and [TT1 the resolution of configuration space sampling 
has been reduced due to the computational cost.) 

Qualitatively, there are no severe changes, although the phase transition line 
widens for higher values of a b . This seems clear, since a wide range of burning prob- 
abilities allows some sub- and supercritical configurations to obtain characteristics 
from the critical case. On the other hand, critical configurations may gain more 
characteristics of either sub- or supercriticality as well, so the transition is smeared 
out over a larger region of parameter space. 

The magnitude of the standard deviation does not change significantly compared 
to simulations performed for a b = (but with otherwise identical parameters). This 
may be somehow surprising, since one could expect the deviation in pbum to result 
in larger differences between simulation runs with the same choice of parameters. 

But of course for the critical or nearly critical case, even with a b = the same 
set of parameters can lead to radically different outcomes. So the deviations can be 
expected to be large already in this case; any change of a b will not significantly in- 
crease them. Nevertheless there are several significant quantitative changes because 
of the modification of the width a b : 

• The burning rates /burn are reduced in the supercritical region, which is clear 
since also in the supercritical case now there will remain several trees with a 
burning probability considerably lower than average. 

• The total burning time tburn is reduced as well as the burning maximum 
^Vburn,max, probably due to similar reasons. 

• While the change in the linear fit parameter ai is not significant, this is not true 
for the exponential fit via a e . Large values of \a e \ indicate a rapid saturation in 
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the cell combination analysis - and such a saturation will definitely be slowed 
down by regions of considerably increased or reduced burning probability. So 
\a e \ is considerably smaller for larger values of ov 
• The dominant Fourier peak is shifted towards lower frequencies while the in- 
tensities remain similar. 

6. Summary and Outlook 

We have studied an enhanced forest fire model that shows a surprisingly sharp 
phase transition in respect to the two main parameters, namely tree density <4tart 
and burning probability pb for a healthy tree with a burning neighbour. 

In a large region of configuration space, the fire (or pathogen) vanishes within 
a few timesteps, whereas in most of the remaining space, the fire affects almost 
the whole population. The border between these two outcomes is a small region of 
"critical" behaviour that can be characterized by cluster formation (on any length 
scale) and exceptional values for many characteristic quantities. (For example the 
total burning time tends to be more than twice as long as in any other case.) 

A redefinition of neighbourhood, controlled by the Moore-parameter ttm did not 
change the qualitative behaviour, but influenced the position of the line of critical 
behaviour in phase space. In contrast, the transition to individual susceptibility had 
only small influence on the global characteristics of the system also on quantitative 
level. 

The model can easily be re-interpreted as a disease spreading model, where sus- 
ceptible, infected and removed individuals replace living, burning and burnt-out 
trees. Our results thus may indicate that the transition from an endemic disease 
(which can only affect a small fraction of the population) to an epidemic (or, de- 
pending on the system, even a pandemic) can be driven by a change of just a few 
percent in one crucial parameter. 

While some analysis is to be done also for this model (like the check whether 
for the critical case a power law arises for the cluster size distribution), there are 
several possible extensions which remain to be investigated. Allowing empty cells 
to be (re)populated changes the model significantly and reduces the importance of 
the initial density. Some results can be found in [2] ; a comprehensive article on this 
topic is in preparation. 
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Fig. 12. Influence of <T(, = 0.1, simulations on a 128 X 128-grid with ttm = 0.5, Apj, = Ad s tart = 
0.025 and 10 iterations: (a) data, (b) deviations. The quantities displayed are the same as in 
figure 1101 
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burning trees: 2nd bin e [0, 1] 



burning trees: 3rd bin e [u, 0.51 667] 







Fig. 13. Influence of <T(, = 0.5, simulations on a 128 X 128-grid with ttm = 0.5, Apj, = Ad s tart = 
0.025 and 10 iterations: (a) data, (b) deviations. The quantities displayed are the same as in 
figure 1101 



